// snow_qhc 合成ims气候场文件
// 像素值为10的时候是积雪


#include <string>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <fstream>
//#include "ogr_spatialref.h"//17
#include "gdal_priv.h"
#include <ctime>
#include <fstream>
#include "../sharedcodes-87/wftools.h"
#include "../sharedcodes-87/ajson5.h"
#include "../sharedcodes-87/wgt.h"


using namespace std;
using namespace ArduinoJson;
using std::cout;


//多期合成分类	IMS原始值
//范围外值	-9999
//缺测	-999
//常年陆地 -900 这个值仅在距平产品中有定义，用于区别距平为0的地区
//一次雪都没有且出现水的情况多于陆地和冰	-907
//一次雪都没有且出现陆地的情况多于水和冰	0（陆地颜色）
//一次雪都没有且出现冰的情况多于水和陆地	-909
//出现一次积雪及以上的像素	>0
#define MDAYS_PIXVAL_OUTSIDE -9999
#define MDAYS_PIXVAL_NODATA -999
#define MDAYS_PIXVAL_LAND_JP   -900
#define MDAYS_PIXVAL_WATER -907
#define MDAYS_PIXVAL_ICE -909


struct CombData{
	void resize(int asize) ;
	vector<int> outspacecount ;
	vector<int> nodatacount ;
	vector<int> snowcount ;
	vector<int> watercount ;
	vector<int> icecount ;
	vector<int> landcount ;
	int goodfile ;
};
void CombData::resize(int asize){
	goodfile = 0 ;
	outspacecount.resize(asize) ;
	nodatacount.resize(asize) ;
	snowcount.resize(asize) ;
	watercount.resize(asize) ;
	icecount.resize(asize) ;
	landcount.resize(asize) ;
	for(int it =0; it <asize;++it ){
		outspacecount[it]=0;
		nodatacount[it]=0;
		snowcount[it] = 0 ;
		watercount[it] = 0 ; 
		icecount[it] = 0 ;
		landcount[it] = 0 ;
	}
}


//计算一组文件的气候场
void computeqhc( vector<string> pickfiles , string outfile )
{
	cout<<"================================"<<endl ;
	CombData cd ;

	bool bProjDataOk = false ;
	double trans[6] ;
	char proj[2048] ; 
	int imgsize = 0 ;
	int imgwid = 0 ;
	int imghei = 0 ;
	int numbadfiles = 0 ;

	wft_descSortStrVec( pickfiles) ;

	for (int i = 0; i < pickfiles.size() ; ++i )
	{//逐个文件处理
		string filepath = pickfiles[i];
		GDALDataset* inds1 = (GDALDataset*) GDALOpen(filepath.c_str() , GA_ReadOnly ) ;
		if( inds1 != 0 ){

			imgwid = inds1->GetRasterXSize() ;
			imghei = inds1->GetRasterYSize() ;
			cout<<filepath << " imgsize "<<imgwid<<","<<imghei<<endl ;
			if( bProjDataOk==false ){
				bProjDataOk = true ;
				inds1->GetGeoTransform(trans);
				strcpy( proj ,  inds1->GetProjectionRef() ) ;
				imgsize = imgwid * imghei ;
				cd.resize(imgsize) ;
			}
			cd.goodfile ++ ;
			//临时缓存
			vector<short> tempindata(imgsize) ;
			//读取当前文件的数据
			inds1->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, imgwid, imghei,
				 tempindata.data() , imgwid, imghei, GDT_Int16, 0, 0, 0);
			for(int it = 0 ; it < imgsize ; ++ it )
			{
				int iv = tempindata[it] ;
				if( iv ==  MDAYS_PIXVAL_OUTSIDE ) {
					cd.outspacecount[it] ++ ;
				}else if( iv == MDAYS_PIXVAL_NODATA ){
					cd.nodatacount[it] ++ ;
				}else if( iv > 0 ){
					cd.snowcount[it] += iv ;
				}else if( iv==0 ){
					cd.landcount[it] ++ ;
				}else if( iv == MDAYS_PIXVAL_ICE ){
					cd.icecount[it]++ ;
				}else if( iv == MDAYS_PIXVAL_WATER ){
					cd.watercount[it] ++ ;
				}
			}
			GDALClose(inds1) ;
		}else{
			numbadfiles ++ ;
		}
	}
	cout<<"100 . "<< endl ;
	cout<<"good file :"<<cd.goodfile<<"; bad file:"<<numbadfiles<<endl ;

	cout<<"saving month qhc ..."<<endl ;

	if(bProjDataOk==false){
		cout<<"Error : geotrans and proj are bad"<<endl ;
		return ;
	}

	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	//temp buffer
	vector<short> tempbuffer( imgsize , 0 ) ;
	//compute
	for(int it = 0 ; it < imgsize; ++ it )//here here 
	{
		if( cd.outspacecount[it] > 0 ){
			tempbuffer[it] = MDAYS_PIXVAL_OUTSIDE ;
		}else if( cd.nodatacount[it] >=4 ){
			tempbuffer[it] = MDAYS_PIXVAL_NODATA ;
		}else if( cd.snowcount[it] > 0 ){
			tempbuffer[it] = cd.snowcount[it]*1.0/ cd.goodfile + 0.5  ;//平均积雪天数
		}else if( cd.watercount[it] > cd.icecount[it] && cd.watercount[it] > cd.landcount[it] ){
			tempbuffer[it] = MDAYS_PIXVAL_WATER ;
		}else if( cd.icecount[it] > cd.watercount[it] && cd.icecount[it] > cd.landcount[it] ){
			tempbuffer[it] = MDAYS_PIXVAL_ICE ;
		}else{
			tempbuffer[it] = 0 ;
		}
	}


	GDALDataset* outds = driver->Create( outfile.c_str() , imgwid,imghei, 1, GDT_Int16, 0);
	outds->SetGeoTransform(trans);
	outds->SetProjection(proj);
	outds->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, imgwid, imghei,
		 tempbuffer.data() , imgwid, imghei, GDT_Int16, 0, 0, 0);			
	GDALClose(outds) ;
	cout<<outfile << " done"<<endl ;
	cout<<"================================"<<endl ;
}


//合成每月气候场 12月 输入文件名 ims_20040201.mon.polar.tif
void ims_month_qhc(const int startyear,const int endyear, string indir )
{
	string inprefix = "ims_";
	string intail = ".mon.polar.tif";
	vector<string> allfiles;

	bool dirok = wft_direxist( indir ) ;
	if( dirok )
	{
		wft_get_allSelectedFiles( indir , inprefix, intail, -1, "", allfiles);
	}else{
		cout<<"Error : indir is not exists."<<indir<<endl ;
		return ;
	}

	vector<string> monthfiles[12] ;
	int n = allfiles.size();
	for (int i = 0; i < n; ++i )
	{//逐个文件处理
		string filepath = allfiles[i];
		string filename = wft_base_name(filepath);
		int year1 = atof( filename.substr(4,4).c_str() ) ; 
		int mon1 = atof( filename.substr(8,2).c_str() ) ; 
		int day1 = atof( filename.substr(10,2).c_str() ) ; 
		//该文件不再计算年份范围内的跳过
		if( year1 < startyear || year1 > endyear ) {
			cout<<"out of year ranges, skip : "<< filename <<endl ;
			continue ;//skip
		}
		int imonth = mon1-1 ;
		monthfiles[imonth].push_back( filepath ) ;
	}
	

	for(int imon = 0 ; imon < 12 ; ++ imon )
	{
		 
		//save IMS_IMS_QHCMON_YYYYMMDD_HHMMSS_VAAAA-BBBB.tif"
		char qhcfile[128] ;
		sprintf(qhcfile,"IMS_IMS_QHCMON_0000%02d01_000000_V%d-%d.tif", imon+1 , startyear , endyear ) ;
		cout<<"begin "<<qhcfile<<endl ;
		computeqhc( monthfiles[imon] , string(qhcfile) ) ;
	}
	cout<<"all month done"<<endl ;
}



//合成每季节气候场   输入文件名 ims_20040301.sea.polar.tif  ims_20050101.yea.polar.tif
void ims_sea_qhc(const int startyear,const int endyear, string indir )
{
	string inprefix = "ims_";
	string intail = ".sea.polar.tif";
	vector<string> allfiles;

	bool dirok = wft_direxist( indir ) ;
	if( dirok )
	{
		wft_get_allSelectedFiles( indir , inprefix, intail, -1, "", allfiles);
	}else{
		cout<<"Error : indir is not exists."<<indir<<endl ;
		return ;
	}
	vector<string> seasonfiles[4] ;
	int n = allfiles.size();
	for (int i = 0; i < n; ++i )
	{//逐个文件处理
		string filepath = allfiles[i];
		string filename = wft_base_name(filepath);
		int year1 = atof( filename.substr(4,4).c_str() ) ; 
		int mon1 = atof( filename.substr(8,2).c_str() ) ; 
		int day1 = atof( filename.substr(10,2).c_str() ) ; 
		//该文件不再计算年份范围内的跳过
		if( year1 < startyear || year1 > endyear ) {
			cout<<"out of year ranges, skip : "<< filename <<endl ;
			continue ;//skip
		}
		int iseason = (mon1-1)/3 ;//3,6,9,12 -> 0,1,2,3
		seasonfiles[iseason].push_back( filepath ) ;
	}
	
	for(int is = 0 ; is < 4 ; ++ is )
	{
		 
		//save IMS_IMS_QHCMON_YYYYMMDD_HHMMSS_VAAAA-BBBB.tif"
		char qhcfile[128] ;
		sprintf(qhcfile,"IMS_IMS_QHCSEA_0000%02d01_000000_V%d-%d.tif", (is+1)*3 , startyear , endyear ) ;
		cout<<"begin "<<qhcfile<<endl ;
		computeqhc( seasonfiles[is] , string(qhcfile) ) ;
	}
	cout<<"all season done"<<endl ;
}



//合成每年气候场 ims_20050101.yea.polar.tif
void ims_yea_qhc( const int startyear,const int endyear,string indir )
{
	string inprefix = "ims_";
	string intail = ".yea.polar.tif";
	vector<string> allfiles;
	bool dirok = wft_direxist( indir ) ;
	if( dirok )
	{
		wft_get_allSelectedFiles( indir , inprefix, intail, -1, "", allfiles);
	}else{
		cout<<"Error : indir is not exists."<<indir<<endl ;
		return ;
	}
	vector<string> pickfiles;
	int n = allfiles.size();
	for (int i = 0; i < n; ++i )
	{//逐个文件处理
		string filepath = allfiles[i];
		string filename = wft_base_name(filepath);
		int year1 = atof( filename.substr(4,4).c_str() ) ; 
		int mon1 = atof( filename.substr(8,2).c_str() ) ; 
		int day1 = atof( filename.substr(10,2).c_str() ) ; 
		//该文件不再计算年份范围内的跳过
		if( year1 < startyear || year1 > endyear ) {
			cout<<"out of year ranges, skip : "<< filename <<endl ;
			continue ;//skip
		}
		pickfiles.push_back( filepath ) ;
	}
	
	//save IMS_IMS_QHCMON_YYYYMMDD_HHMMSS_VAAAA-BBBB.tif"
	char qhcfile[128] ;
	sprintf(qhcfile,"IMS_IMS_QHCYEA_00000101_000000_V%d-%d.tif",  startyear , endyear ) ;
	cout<<"begin "<<qhcfile<<endl ;
	computeqhc( pickfiles, string(qhcfile) ) ;

	cout<<"all season done"<<endl ;
}






int main(int argc, char** argv)
{
	std::cout << "compute ims qhc."<<endl ;
	cout<<"usage:ims_qhc config.json"<<endl ;
	cout<<"config.json:"<<endl ;
	cout<<"  {\"startyear\":2005,\"endyear\":2020,"<<endl ;
	cout<<"   \"mondir\":\"/somedir\","<<endl ;
	cout<<"   \"seadir\":\"/somedir\","<<endl ;
	cout<<"   \"yeadir\":\"/somedir\" "<<endl ;
	cout<<"   }"<<endl ;
	cout<<"input filename should be ims_20040201.mon.polar.tif"<<endl ;
	cout<<"output filename should be IMS_IMS_QHCMON_YYYYMMDD_HHMMSS_VAAAA-BBBB.tif"<<endl ;
	cout<<"AAAA start year, BBBB end year."<<endl ;
	cout<<"v2.0 2020-12-5"<<endl ;
	cout<<"v2.1 use config.json 2021-3-9"<<endl ;
	cout<<"v2.2 use json dir, update method. 2021-3-17"<<endl ;
	cout<<"v2.2.1 sort in files 2021-3-18"<<endl ;

	if( argc != 2 ){
		cout<<"params not 2"<<endl ;
		return 11 ;
	}
	//
	cout<<"gdal register"<<endl ;
	GDALAllRegister();

	//
	string configfile = argv[1] ;
	cout<<"parseing "<<configfile<<endl ;
	
	std::ifstream ifs(configfile.c_str());
	DynamicJsonBuffer jsonBuffer;
	JsonObject& root = jsonBuffer.parseObject(ifs);
	int startYear = root["startyear"].as<int>() ;
	int endYear = root["endyear"].as<int>() ;
	string mondir = root["mondir"].as<char*>() ;
	string seadir = root["seadir"].as<char*>() ;
	string yeadir = root["yeadir"].as<char*>() ;

    ims_month_qhc( startYear , endYear  , mondir   ) ;

    ims_sea_qhc( startYear , endYear, seadir ) ;

    ims_yea_qhc( startYear , endYear , yeadir  ) ;


	std::cout << "All done." << std::endl;

	return 0;
}

